function Calibration_SZHW
clc;clf;close all;

% Created on Thu Jan 03 2019

% The SZHW model and implied volatilities; The calibration takes place for 
% a single maturity time. The numbers reported in the table may differ with respect to 
% function outputs. Here we assume that the market option prices are given by the
% Heston model while in the book we consider real market data.

CP  = 'c';  
        
% HW model parameter settings

lambda = 1.12;
eta   = 0.02;
S0    = 100.0;
T     = 10.0;     
r     = 0.033;
   
% Fixed SZHW model parameters

kappa = 0.5;
Rxr   = 0.;
    
% Strike range

K = linspace(50.0,200.0,20)';

% We define a ZCB curve (obtained from the market)

P0T = @(T) exp(-r*T);
       
% Settings for the COS method

N = 2000;
L = 15;  

% ################## Here we define market call option prices #################
% We assume, in this experiment that the reference prices are generated by the Heston model
% Note that for the Heston model we deal with the discounted ChF, i.e. we have an
% additional -r*tau term in the exponent

% ChF for the Heston model

kappaH = 0.5;
vbarH  = 0.0770;
v0H    = 0.0107;
rhoH   = -0.6622;
gammaH = 0.35;
cfHeston = @(u)ChFHestonModelDisc(u, T, kappaH,vbarH,gammaH,rhoH, v0H, r);
referencePrice = CallPutOptionPriceCOSMthd_StochIR(cfHeston, CP, S0, T, K, N, L,P0T(T));

figure(1); hold on;
plot(K,referencePrice)
xlabel('strike,K')
ylabel('Optio price')
grid on

% Global search algorithm
% The search will terminate after running 30 seconds 

gs = GlobalSearch('FunctionTolerance',1e-3,'NumTrialPoints',300,'XTolerance',1e-3,'MaxTime',30);

% Define a target function

trgt        = @(x)TargetVal(CP,kappa,x(1),x(2),Rxr,x(3),x(4),x(5),eta,lambda,K,referencePrice,S0,T,P0T);
lowerBound  = [0.1, 0.01,-0.85, -0.85, 0.001];
upperBound  = [0.8,  0.6, 0.85, -0.0, 0.8];
x0          = [0.1, 0.3, 0.0,-0.5,0.3];
problemSZHW = createOptimProblem('fmincon','x0',x0,'objective',trgt,'lb',lowerBound,'ub',upperBound);
x           = run(gs,problemSZHW);

% Optimal set of parameters

gamma_opt = x(1);
sigmaBar_opt= x(2);
Rrsigma_opt = x(3);
Rxsigma_opt = x(4);
sigma0_opt = x(5);

% Plot final result

cf = @(u)ChFSZHW(u,P0T,T,kappa,sigmaBar_opt,gamma_opt,lambda,eta,Rxsigma_opt,Rxr,Rrsigma_opt,sigma0_opt);
valCOS    = CallPutOptionPriceCOSMthd_StochIR(cf,CP,S0,T,K,N,L,P0T(T));
plot(K,valCOS,'--r')
legend('reference option prices','calibrated option prices')
xlabel('strike')
ylabel('option price')

fprintf('Optimal parameters are: gamma = %.4f, sigmaBar = %.4f, Rrsigma = %.4f, Rxsigma = %.4f, sigma0 = %.4f ', x(1),x(2),x(3),x(4),x(5))

function value = TargetVal(CP,kappa,gamma,sigmaBar,Rxr,Rrsigma,Rxsigma,sigma0,eta,lambda,K,marketPrice,S0,T,P0T)

% Settings for the COS method

N = 1000;
L = 15;
cf = @(u)ChFSZHW(u,P0T,T,kappa, sigmaBar,gamma,lambda,eta,Rxsigma,Rxr,Rrsigma,sigma0);
valCOS    = CallPutOptionPriceCOSMthd_StochIR(cf,CP,S0,T,K,N,L,P0T(T));

% Error is a squared sum

value = sum((marketPrice-valCOS).^2);
fprintf('Total Error = %.4f',value)
fprintf('\n')

function cf=ChFHestonModelDisc(u, tau, kappa,vBar,gamma,rho, v0, r)
i     = complex(0,1);

% Functions D_1 and g

D_1  = sqrt(((kappa -i*rho*gamma.*u).^2+(u.^2+i*u)*gamma^2));
g    = (kappa- i*rho*gamma*u-D_1)./(kappa-i*rho*gamma*u+D_1);    

% Complex-valued functions A and C

C = (1/gamma^2)*(1-exp(-D_1*tau))./(1-g.*exp(-D_1*tau)).*(kappa-gamma*rho*i*u-D_1);
A = -r*tau + i*u*r*tau + kappa*vBar*tau/gamma^2 * (kappa-gamma*rho*i*u-D_1)-2*kappa*vBar/gamma^2*log((1-g.*exp(-D_1*tau))./(1-g));

% ChF for the Heston model

cf = exp(A + C * v0);

function value = ChFSZHW(u,P0T,tau,kappa, sigmabar,gamma,lambda,eta,Rxsigma,Rxr,Rrsigma,sigma0)
v_D = D(u,tau,kappa,Rxsigma,gamma);
v_E = E(u,tau,lambda,gamma,Rxsigma,Rrsigma,Rxr,eta,kappa,sigmabar);
v_A = A(u,tau,P0T,eta,lambda,Rxsigma,Rrsigma,Rxr,gamma,kappa,sigmabar);

% Initial variance 

v0     =sigma0^2;

% Characteristic function of the SZHW model

value = exp(v0*v_D + sigma0*v_E + v_A);
        
function value=C(u,tau,lambda)
    i     = complex(0,1);
    value = 1/lambda*(i*u-1).*(1-exp(-lambda*tau)); 

function value=D(u,tau,kappa,Rxsigma,gamma)
    i=complex(0,1);
    a_0=-1/2*u.*(i+u);
    a_1=2*(gamma*Rxsigma*i*u-kappa);
    a_2=2*gamma^2;
    d=sqrt(a_1.^2-4*a_0.*a_2);
    g=(-a_1-d)./(-a_1+d);    
value=(-a_1-d)./(2*a_2*(1-g.*exp(-d*tau))).*(1-exp(-d*tau));

function value=E(u,tau,lambda,gamma,Rxsigma,Rrsigma,Rxr,eta,kappa,sigmabar)
    i=complex(0,1);
    a_0=-1/2*u.*(i+u);
    a_1=2*(gamma*Rxsigma*i*u-kappa);
    a_2=2*gamma^2;
    d  =sqrt(a_1.^2-4*a_0.*a_2);
    g =(-a_1-d)./(-a_1+d);    
    
    c_1=gamma*Rxsigma*i*u-kappa-1/2*(a_1+d);
    f_1=1./c_1.*(1-exp(-c_1*tau))+1./(c_1+d).*(exp(-(c_1+d)*tau)-1);
    f_2=1./c_1.*(1-exp(-c_1*tau))+1./(c_1+lambda).*(exp(-(c_1+lambda)*tau)-1);
    f_3=(exp(-(c_1+d)*tau)-1)./(c_1+d)+(1-exp(-(c_1+d+lambda)*tau))./(c_1+d+lambda);
    f_4=1./c_1-1./(c_1+d)-1./(c_1+lambda)+1./(c_1+d+lambda);
    f_5=exp(-(c_1+d+lambda).*tau).*(exp(lambda*tau).*(1./(c_1+d)-exp(d*tau)./c_1)+exp(d*tau)./(c_1+lambda)-1./(c_1+d+lambda)); 

    I_1=kappa*sigmabar./a_2.*(-a_1-d).*f_1;
    I_2=eta*Rxr*i*u.*(i*u-1)./lambda.*(f_2+g.*f_3);
    I_3=-Rrsigma*eta*gamma./(lambda.*a_2).*(a_1+d).*(i*u-1).*(f_4+f_5);
    value=exp(c_1*tau).*1./(1-g.*exp(-d*tau)).*(I_1+I_2+I_3);
    
function value=A(u,tau,P0T,eta,lambda,Rxsigma,Rrsigma,Rxr,gamma,kappa,sigmabar)
    i=complex(0,1);
    a_0=-1/2*u.*(i+u);
    a_1=2*(gamma*Rxsigma*i*u-kappa);
    a_2=2*gamma^2;
    d  =sqrt(a_1.^2-4*a_0.*a_2);
    g =(-a_1-d)./(-a_1+d); 
    f_6=eta^2/(4*lambda^3)*(i+u).^2*(3+exp(-2*lambda*tau)-4*exp(-lambda*tau)-2*lambda*tau);
    A_1=1/4*((-a_1-d)*tau-2*log((1-g.*exp(-d*tau))./(1-g)))+f_6;
   
    % Integration within the function A(u,tau)

    value=zeros(1,length(u));   
    N=500;
    arg=linspace(0,tau,N);

    % Solve the integral for A

    for k=1:length(u)
       E_val=E(u(k),arg,lambda,gamma,Rxsigma,Rrsigma,Rxr,eta,kappa,sigmabar);
       C_val=C(u(k),arg,lambda);
       f=(kappa*sigmabar+1/2*gamma^2*E_val+gamma*eta*Rrsigma*C_val).*E_val;
       value(1,k)=trapz(arg,f);
    end
    value = value + A_1;   
    
% Correction for the interest rate term structure    

help = eta^2/(2*lambda^2)*(tau+2/lambda*(exp(-lambda*tau)-1)-1/(2*lambda)*(exp(-2*lambda*tau)-1));
correction = (i*u-1).*(log(1/P0T(tau))+help);
value = value + correction;

function value = CallPutOptionPriceCOSMthd_StochIR(cf,CP,S0,tau,K,N,L,P0T)
i = complex(0,1);

% cf   - Characteristic function, in the book denoted as \varphi
% CP   - C for call and P for put
% S0   - Initial stock price
% r    - Interest rate (constant)
% tau  - Time to maturity
% K    - Vector of strike prices
% N    - Number of expansion terms
% L    - Size of truncation domain (typ.:L=8 or L=10)

x0 = log(S0 ./ K);   

% Truncation domain

a = 0 - L * sqrt(tau); 
b = 0 + L * sqrt(tau);

k = 0:N-1;             % Row vector, index for expansion terms
u = k * pi / (b - a);  % ChF arguments

H_k = CallPutCoefficients('P',a,b,k);
temp = (cf(u) .* H_k).';
temp(1) = 0.5 * temp(1);      % Multiply the first element by 1/2

mat = exp(i * (x0 - a) * u);  % Matrix-vector manipulations

% Final output

value = K .* real(mat * temp);

% Use the put-call parity to determine call prices (if needed)

if lower(CP) == 'c' || CP == 1
 value = value + S0 - K * P0T; 
end

% Coefficients H_k for the COS method

function H_k = CallPutCoefficients(CP,a,b,k)
 if lower(CP) == 'c' || CP == 1
  c = 0;
  d = b;
  [Chi_k,Psi_k] = Chi_Psi(a,b,c,d,k);
   if a < b && b < 0.0
   H_k = zeros([length(k),1]);
   else
   H_k = 2.0 / (b - a) * (Chi_k - Psi_k);
   end
 elseif lower(CP) == 'p' || CP == -1
  c = a;
  d = 0.0;
  [Chi_k,Psi_k]  = Chi_Psi(a,b,c,d,k);
   H_k = 2.0 / (b - a) * (- Chi_k + Psi_k);    
 end

function [chi_k,psi_k] = Chi_Psi(a,b,c,d,k)
 psi_k  = sin(k * pi * (d - a) / (b - a)) - sin(k * pi * (c - a)/(b - a));
 psi_k(2:end) = psi_k(2:end) * (b - a) ./ (k(2:end) * pi);
 psi_k(1)  = d - c;
 
 chi_k = 1.0 ./ (1.0 + (k * pi / (b - a)).^2); 
 expr1 = cos(k * pi * (d - a)/(b - a)) * exp(d)  - cos(k * pi... 
      * (c - a) / (b - a)) * exp(c);
 expr2 = k * pi / (b - a) .* sin(k * pi * ...
      (d - a) / (b - a))   - k * pi / (b - a) .* sin(k... 
      * pi * (c - a) / (b - a)) * exp(c);
 chi_k = chi_k .* (expr1 + expr2);
 
